#include "AnaDef.h"
#include "AnaAliDef.h"
#include "XGLUtils.h"

void MIV0()
{
  TFile f("/hera/alice/miranov/highPt/LHC13c/pass2/V0Selected.root");

  const TString tn[]={"K0","Lambda","ALambda","Gamma"};
  const Double_t mass[]={0.139570, 0.938272,0.938272, 5.109989e-04};
  const Int_t ntree = sizeof(tn)/sizeof(TString);
  TList *ll=new TList;

  for(Int_t itree=0; itree<ntree; itree++){
    TTree * tt = (TTree*) gDirectory->Get(Form("tree%s",tn[itree].Data()));
    if(!tt){
      printf("tt null!! %s\n", tn[itree].Data()); 
      gDirectory->ls();
      exit(1);  
    }

    for(Int_t itk=0; itk<2; itk++){
      if(1){
        const TString hn = Form("h%strack%d", tn[itree].Data(), itk);
        TH2D *hh = new TH2D(hn,hn,400, 0.1, 10, 400, 0, 10);
        XGLUtils::BinLogX(hh->GetXaxis());
        ll->Add(hh);
        
        const TString var = Form("track%d.fTRDsignal:track%d.P()", itk, itk);
        const TString cut = Form("track%d.fTRDNchamberdEdx==6 && track%d.fTRDNclusterdEdx/track%d.fTRDNchamberdEdx>17", itk, itk, itk);
        hh->SetTitle((var+" ")+cut);
        
        tt->Draw((var+">>")+hn, cut);
      }
      if(1){
        const TString hn = Form("bg%strack%d", tn[itree].Data(), itk);
        TH2D *hh = new TH2D(hn,hn,500, 0.1, 1e4, 400, 0, 10);
        XGLUtils::BinLogX(hh->GetXaxis());
        ll->Add(hh);
        
        const TString var = Form("track%d.fTRDsignal:track%d.P()/%e", itk, itk, mass[itree]);
        const TString cut = Form("track%d.fTRDNchamberdEdx==6 && track%d.fTRDNclusterdEdx/track%d.fTRDNchamberdEdx>17", itk, itk, itk);
        hh->SetTitle((var+" ")+cut);
        
        tt->Draw((var+">>")+hn, cut);
      }
    }
  }

  const Double_t thres = 10;

  if(1){
    TH2D *hpion = (TH2D*)((TH2D*) ll->FindObject("hK0track0"))->Clone("hpion");
    hpion->Add((TH2D*) ll->FindObject("hK0track1")); ll->Add(hpion); 
    //hpion->RebinX(4);
    //hpion->RebinY(4);
    hpion->GetXaxis()->SetRangeUser(0.3, 5);
    hpion->GetYaxis()->SetRangeUser(0,6);
    
    TH2D *htot = (TH2D*) hpion->Clone("htot"); ll->Add(htot);
    TH2D *hnorpion=XGLUtils::NormalHist(hpion, thres, kTRUE);ll->Add(hnorpion);
    
    TH2D *hproton = (TH2D*)((TH2D*) ll->FindObject("hLambdatrack0"))->Clone("hproton");
    hproton->Add((TH2D*) ll->FindObject("hALambdatrack1"));ll->Add(hproton); 
    //hproton->RebinX(4);
    //hproton->RebinY(4);
    
    htot->Add(hproton);
    TH2D *hnorproton=XGLUtils::NormalHist(hproton, thres, kTRUE);ll->Add(hnorproton);
    
    TH2D *helectron = (TH2D*)((TH2D*) ll->FindObject("hGammatrack0"))->Clone("helectron");
    helectron->Add((TH2D*) ll->FindObject("hGammatrack1"));ll->Add(helectron); 
    //helectron->RebinX(4);
    //helectron->RebinY(4);
    
    htot->Add(helectron);
    TH2D *hnorelectron=XGLUtils::NormalHist(helectron, thres, kTRUE);ll->Add(hnorelectron);
    
    TH2D *hnortot=XGLUtils::NormalHist(htot, thres, kFALSE);ll->Add(hnortot);
  }
  if(1){
    TH2D *bgpion = (TH2D*)((TH2D*) ll->FindObject("bgK0track0"))->Clone("bgpion");
    bgpion->Add((TH2D*) ll->FindObject("bgK0track1")); ll->Add(bgpion); 
    //bgpion->RebinX(4);
    //bgpion->RebinY(4);
    //bgpion->GetXaxis()->SetRangeUser(0.3, 5);
    bgpion->GetYaxis()->SetRangeUser(0,6);
    
    TH2D *htot = (TH2D*) bgpion->Clone("bgtot"); ll->Add(htot);
    TH2D *hnorpion=XGLUtils::NormalHist(bgpion, thres, kTRUE);ll->Add(hnorpion);
    
    TH2D *bgproton = (TH2D*)((TH2D*) ll->FindObject("bgLambdatrack0"))->Clone("bgproton");
    bgproton->Add((TH2D*) ll->FindObject("bgALambdatrack1"));ll->Add(bgproton); 
    //bgproton->RebinX(4);
    //bgproton->RebinY(4);
    
    //htot->Add(bgproton);
    TH2D *hnorproton=XGLUtils::NormalHist(bgproton, thres, kTRUE);ll->Add(hnorproton);
    
    TH2D *bgelectron = (TH2D*)((TH2D*) ll->FindObject("bgGammatrack0"))->Clone("bgelectron");
    bgelectron->Add((TH2D*) ll->FindObject("bgGammatrack1"));ll->Add(bgelectron); 
    //bgelectron->RebinX(4);
    //bgelectron->RebinY(4);
    
    htot->Add(bgelectron);
    TH2D *hnorelectron=XGLUtils::NormalHist(bgelectron, thres, kTRUE);ll->Add(hnorelectron);
    
    TH2D *hnortot=XGLUtils::NormalHist(htot, thres, kFALSE);ll->Add(hnortot);

    TH1D * hamp = 0x0; TH1D * hmpv = 0x0; TH1D * hwid = 0x0; TH1D * hchi = 0x0; TH1D * hres = 0x0;
    XGLUtils::FitSlicesY(hnortot, /*XGLUtils::RawGaus*/ 0x0, hamp, hmpv, hwid, hres, hchi, 1, 0,1,0x0); 
    ll->Add(hamp);
    ll->Add(hmpv);
    ll->Add(hwid);
    ll->Add(hchi);
    ll->Add(hres);
  }

  TFile *fout=new TFile("outMIV0.root","recreate");
  ll->Write();
  fout->Save();
  fout->Close();
  delete fout;
}

int main()
{
  MIV0();
  return 0;
}
